#-----------------------------------------------------------------------------#
#### IMPORT ####
#-----------------------------------------------------------------------------#
library("lfe")
library(texreg)
library(data.table)
library(readstata13)
# library(plm)
library(Hmisc)
library(haven)
library(plyr,include.only = "rbind.fill")
library(R2HTML)
memory.limit(size=64000)

# Functions
simplelag <- function (x,by=NULL,mylag=1,outside=NA) {
  myend <- length(x)- mylag
  if(!is.null(by))  {
    lby <- c(replicate(mylag,""),as.character(by[1:myend]))
    y0 <- c(replicate(mylag,outside),x[1:myend])
    y <- ifelse(as.character(by)==lby,y0,outside)
  }
  else {
    y <- c(replicate(mylag,outside),x[1:myend])
  }
}



forward <- function (x,by=NULL,myforward=1,outside=NA) {
  mystart <- myforward+1
  if(!is.null(by))  {
    fby <- c(as.character(by[mystart:length(x)],replicate(myforward,"")))
    y0 <- c(x[mystart:length(x)],replicate(myforward,outside))
    y <- ifelse(as.character(by)==fby,y0,outside)
  }
  else {
    y <- c(x[mystart:length(x)],replicate(myforward,outside))
  }
}


retain <- function(x,event,outside=NA)
{
  indices <- c(1,which(event==TRUE),length(x)+1)
  values <- c(outside,x[event %in% TRUE])
  y <- rep(values,diff(indices))
}

obtain <- function(x,event,outside=NA)
{
  indices <- c(0,which(event==TRUE),length(x))
  values <- c(x[event %in% TRUE],outside)
  y <- rep(values,diff(indices))
}

wtd.summary <- function (x,w=NA) {
  name <- deparse(substitute(x))
  nb_obs <- length(x)
  sum_wgt <- sum(w[is.na(x) ==F])
  nb_missing <- sum(is.na(x))
  sum_wgt_missing <- sum(w[is.na(x) ==T])
  wtd_mean <- weighted.mean(x,w=w,na.rm=T)
  wtd_sd <- wtd.var(x,w=w,na.rm=T)**0.5
  min <- min(x,na.rm=T)
  p1 <- wtd.quantile(x,weights=w,probs=0.01)
  p10 <- wtd.quantile(x,weights=w,probs=0.10)
  q1 <- wtd.quantile(x,weights=w,probs=0.25)
  q2 <- wtd.quantile(x,weights=w,probs=0.50)
  q3 <- wtd.quantile(x,weights=w,probs=0.75)
  p90 <- wtd.quantile(x,weights=w,probs=0.90)
  p99 <- wtd.quantile(x,weights=w,probs=0.99)
  max <- max(x,na.rm=T)
  structure(data.frame(name,nb_obs,sum_wgt,nb_missing,sum_wgt_missing,wtd_mean,wtd_sd,min,p1,p10,q1,q2,q3,p90,p99,max))
}

#-----------------------------------------------------------------------------#
#### DATA IMPORT ####
#-----------------------------------------------------------------------------#

setwd("C:\\Users\\Public\\Documents\\Olivier\\Results\\Establishment_Composition")
gg <- data.table(readRDS("gg.rds"))
gg <- gg[order(gg$est,gg$year),]
gg[,c("sf9910","f9910_w","f9910xf9910","f9099","f9910","sf9010","lndnbwkrs_neg",
"lcount","lnbwkrs"):=NULL]

#-----------------------------------------------------------------------------#
#### DATA MANAGEMENT ####
#-----------------------------------------------------------------------------#


gg$f9010_min <- ave(gg$f9010,gg$est,FUN=min)
# table(gg$f9010_min)
field <- gg$f9010_w2>0 & is.na(gg$f9010_w2)==F & gg$f9010_min>1
# table(field)

gg <- gg[field==T,]
rownames(gg) <- NULL

str(gg)
gg$nblayoff_e <- ave(gg$nblayoffs>0,gg$est,FUN=sum)
gg$nblayoff_o <- ave(gg$nblayoffs>0,gg$est,FUN=cumsum)

gg$nbyear_e <- ave(gg$year>0,gg$est,FUN=sum)
gg$nbyear_o <- ave(gg$year>0,gg$est,FUN=cumsum)


gg$playoff <- retain(gg$year*(gg$nblayoffs>0),gg$nblayoffs>0 | gg$est!=simplelag(gg$est),outside=0)
gg$flayoff <- obtain(gg$year*(gg$nblayoffs>0),gg$nblayoffs>0  | gg$est != forward(gg$est),outside=0)

gg$playoff <- ifelse(gg$nblayoffs>0,simplelag(gg$playoff,by=gg$est,outside=0),gg$playoff)
gg$flayoff <- ifelse(gg$nblayoffs>0,forward(gg$flayoff,by=gg$est,outside=0),gg$flayoff)

# gg$playoff_sz <- retain(gg$nblayoffs,gg$nblayoffs>0 | gg$est!=simplelag(gg$est),outside=0)
# gg$flayoff_sz <- obtain(gg$nblayoffs,gg$nblayoffs>0  | gg$est != forward(gg$est),outside=0)


# Stacking data together
gc()
fff <- NULL
for (y in c(2002:2014)) {
  gg$treated_y <- ( gg$year<y+5 &  gg$year>y-5 &
                    ((gg$nblayoffs>0)*(gg$year==y) | gg$playoff==y | gg$flayoff==y) &
                    (gg$flayoff==0 | ( ( (gg$flayoff>y)*gg$year ) < gg$flayoff )) &
                    (gg$playoff==0 | (gg$year > ((gg$playoff<y)*gg$playoff)))  &
                    ((gg$nblayoffs>0)*(gg$year != y) ==0)
                      )*1
  gg$treated_y_sum <- ave(gg$treated_y,gg$est,FUN=sum)
  
  gg$control_y <- ( gg$treated_y_sum==0 &
                      gg$year<y+5 &  gg$year>y-5 &
                      (gg$nblayoffs==0) &
                      (gg$flayoff==0 | (gg$year < gg$flayoff) ) &
                      (gg$playoff==0 | (gg$year > gg$playoff  )) 
                    )*1
  
  ff <- gg[gg$treated_y %in% 1 | gg$control_y %in% 1,
           c("est","year","count","lnnbwkrs_cumneg","f9010xf9010","f9010_w","f9010_w2","nblayoffs",
             "playoff","flayoff","r_lyof","treated_y","control_y")]
  ff$eyear <- y
  ff$y_nblayoffs <- ave((ff$nblayoffs>0)*ff$year,ff$est,FUN=function(x) max(x,na.rm=T))
  ff$y_r_lyof <- ave(ff$r_lyof,ff$est,FUN=function(x) max(x,na.rm=T))
  ff$min_playoff <- ave(ff$playoff,ff$est,FUN=function(x) min(x,na.rm=T))
  ff$nb_obs <- ave(ff$year>0,ff$est,FUN=sum)
  ff$min_year <- ave(ff$year,ff$est,FUN=min)
  ff$max_year <- ave(ff$year,ff$est,FUN=max)
  fff <- rbind(fff,ff[ff$nb_obs>1 &  ff$min_year<=y & ff$max_year>=y,])
}



table(fff$y_nblayoffs,fff$treated_y)
nrow(fff)

table(fff$nblayoffs>0,fff$year)

#-----------------------------------------------------------------------------#
#### DESCRIPTIVES ####
#-----------------------------------------------------------------------------#
tt1 <- table(fff$year[fff$year==fff$eyear],fff$treated_y[fff$year==fff$eyear])
tt2 <- wtd.table(paste(fff$year[fff$year==fff$eyear],fff$treated_y[fff$year==fff$eyear]),weights=fff$f9010_w2[fff$year==fff$eyear])
data.frame(tt2)
tt1
file.remove("layoff_stacked_des.html")
HTML(tt1,"layoff_stacked_des.html")
HTML(tt2,"layoff_stacked_des.html")

#-----------------------------------------------------------------------------#
#### MODELS ####
#-----------------------------------------------------------------------------#
fff$firm <- substr(fff$est,1,9)
gc()
mod90_3 <- felm(I(100*f9010xf9010) ~ factor(paste(eyear,year))
                + I((year - y_nblayoffs) %in% -4 )
                + I((year - y_nblayoffs) %in% -3 )
                + I((year - y_nblayoffs) %in% -2 )
                # + I((year - y_nblayoffs) %in% -1 & y_nblayoffs==eyear)
                + I((year - y_nblayoffs) %in% 0 )
                + I((year - y_nblayoffs) %in% 1 )
                + I((year - y_nblayoffs) %in% 2 )
                + I((year - y_nblayoffs) %in% 3 )
                + I((year - y_nblayoffs) %in% 4 )
                |paste(eyear,est)|0|firm, 
                weights=fff$f9010_w2, 
                data=fff  )
summary(mod90_3)
screenreg(list(mod90_3),stars=c(0.1,0.05,0.01),digits=3)
htmlreg(list(mod90_3),stars=c(0.1,0.05,0.01),file="did_layoff_stacked.html",digits=3)



mod90_3 <- felm(I(100*f9010xf9010) ~ factor(paste(eyear,year))
                + I((year - y_nblayoffs) %in% -4 )
                + I((year - y_nblayoffs) %in% -3 )
                + I((year - y_nblayoffs) %in% -2 )
                # + I((year - y_nblayoffs) %in% -1 & y_nblayoffs==eyear)
                + I((year - y_nblayoffs) %in% 0 )
                + I((year - y_nblayoffs) %in% 1 )
                + I((year - y_nblayoffs) %in% 2 )
                + I((year - y_nblayoffs) %in% 3 )
                + I((year - y_nblayoffs) %in% 4 )
                + log(count)
                + lnnbwkrs_cumneg
                |paste(eyear,est)|0|firm, 
                weights=fff$f9010_w2, 
                data=fff  )
summary(mod90_3)
screenreg(list(mod90_3),stars=c(0.1,0.05,0.01),digits=3)
htmlreg(list(mod90_3),stars=c(0.1,0.05,0.01),file="did_layoff_stacked_size.html",digits=3)

gc()
mod90_3 <- felm(I(100*f9010xf9010) ~ factor(paste(eyear,year))
                + I(y_r_lyof*((year - y_nblayoffs) %in% -4 ))
                + I(y_r_lyof*((year - y_nblayoffs) %in% -3 ))
                + I(y_r_lyof*((year - y_nblayoffs) %in% -2 ))
                # + I(y_r_lyof*((year - y_nblayoffs) %in% -1 & y_nblayoffs==eyear))
                + I(y_r_lyof*((year - y_nblayoffs) %in% 0 ))
                + I(y_r_lyof*((year - y_nblayoffs) %in% 1 ))
                + I(y_r_lyof*((year - y_nblayoffs) %in% 2 ))
                + I(y_r_lyof*((year - y_nblayoffs) %in% 3 ))
                + I(y_r_lyof*((year - y_nblayoffs) %in% 4 ))
                |paste(eyear,est)|0|firm, 
                weights=fff$f9010_w2, 
                data=fff  )
summary(mod90_3)
screenreg(list(mod90_3),stars=c(0.1,0.05,0.01),digits=3)
htmlreg(list(mod90_3),stars=c(0.1,0.05,0.01),file="did_r_lyof_stacked.html",digits=3)

mod90_3 <- felm(I(100*f9010xf9010) ~ factor(paste(eyear,year))
                + I(y_r_lyof*((year - y_nblayoffs) %in% -4 ))
                + I(y_r_lyof*((year - y_nblayoffs) %in% -3 ))
                + I(y_r_lyof*((year - y_nblayoffs) %in% -2 ))
                # + I(y_r_lyof*((year - y_nblayoffs) %in% -1 & y_nblayoffs==eyear))
                + I(y_r_lyof*((year - y_nblayoffs) %in% 0 ))
                + I(y_r_lyof*((year - y_nblayoffs) %in% 1 ))
                + I(y_r_lyof*((year - y_nblayoffs) %in% 2 ))
                + I(y_r_lyof*((year - y_nblayoffs) %in% 3 ))
                + I(y_r_lyof*((year - y_nblayoffs) %in% 4 ))
                + log(count)
                + lnnbwkrs_cumneg
                |paste(eyear,est)|0|firm, 
                weights=fff$f9010_w2, 
                data=fff  )
summary(mod90_3)
screenreg(list(mod90_3),stars=c(0.1,0.05,0.01),digits=3)
htmlreg(list(mod90_3),stars=c(0.1,0.05,0.01),file="did_r_lyof_stacked_size.html",digits=3)
gc()


#------------------------------------#
#### Trend analysis  ####
#------------------------------------#
gc()
mod90_7 <- felm(I(100*f9010xf9010) ~ year
                # +factor(eyear):year
                # + I((year - y_nblayoffs) %in% -4 )
                # + I((year - y_nblayoffs) %in% -3 )
                # + I((year - y_nblayoffs) %in% -2 )
                # # + I((year - y_nblayoffs) %in% -1 & y_nblayoffs==eyear)
                # + I((year - y_nblayoffs) %in% 0 )
                # + I((year - y_nblayoffs) %in% 1 )
                # + I((year - y_nblayoffs) %in% 2 )
                # + I((year - y_nblayoffs) %in% 3 )
                # + I((year - y_nblayoffs) %in% 4 )
                |paste(eyear,est)|0|firm, 
                weights=fff$f9010_w2, 
                data=fff  )
summary(mod90_7)

screenreg(list(mod90_7),stars=c(0.1,0.05,0.01),digits=5)

mod90_8 <- felm(I(100*f9010xf9010) ~ year
                # +factor(eyear):year
                + I((year - y_nblayoffs) %in% -4 )
                + I((year - y_nblayoffs) %in% -3 )
                + I((year - y_nblayoffs) %in% -2 )
                # + I((year - y_nblayoffs) %in% -1 & y_nblayoffs==eyear)
                + I((year - y_nblayoffs) %in% 0 )
                + I((year - y_nblayoffs) %in% 1 )
                + I((year - y_nblayoffs) %in% 2 )
                + I((year - y_nblayoffs) %in% 3 )
                + I((year - y_nblayoffs) %in% 4 )
                |paste(eyear,est)|0|firm, 
                weights=fff$f9010_w2, 
                data=fff  )
summary(mod90_8)
screenreg(list(mod90_8),stars=c(0.1,0.05,0.01),digits=5)


mod90_9 <- felm(I(100*f9010xf9010) ~ year
                + I(y_r_lyof*((year - y_nblayoffs) %in% -4 ))
                + I(y_r_lyof*((year - y_nblayoffs) %in% -3 ))
                + I(y_r_lyof*((year - y_nblayoffs) %in% -2 ))
                # + I(y_r_lyof*((year - y_nblayoffs) %in% -1 & y_nblayoffs==eyear))
                + I(y_r_lyof*((year - y_nblayoffs) %in% 0 ))
                + I(y_r_lyof*((year - y_nblayoffs) %in% 1 ))
                + I(y_r_lyof*((year - y_nblayoffs) %in% 2 ))
                + I(y_r_lyof*((year - y_nblayoffs) %in% 3 ))
                + I(y_r_lyof*((year - y_nblayoffs) %in% 4 ))
                |paste(eyear,est)|0|firm, 
                weights=fff$f9010_w2, 
                data=fff  )
summary(mod90_9)
screenreg(list(mod90_7,mod90_8,mod90_9),stars=c(0.1,0.05,0.01),digits=5)
htmlreg(list(mod90_7,mod90_8,mod90_9),stars=c(0.1,0.05,0.01),file="did_lyof_stacked_trend.html",digits=5)



#-----------------------------------------------------------------------------#
#### DATA Layoffs>5 ####
#-----------------------------------------------------------------------------#
hh <- gg

str(hh)
hh$nblayoff_e <- ave(hh$nblayoffs>5,hh$est,FUN=sum)
hh$nblayoff_o <- ave(hh$nblayoffs>5,hh$est,FUN=cumsum)

hh$nbyear_e <- ave(hh$year>0,hh$est,FUN=sum)
hh$nbyear_o <- ave(hh$year>0,hh$est,FUN=cumsum)


hh$playoff <- retain(hh$year*(hh$nblayoffs>5),hh$nblayoffs>5 | hh$est!=simplelag(hh$est),outside=0)
hh$flayoff <- obtain(hh$year*(hh$nblayoffs>5),hh$nblayoffs>5  | hh$est != forward(hh$est),outside=0)

hh$playoff <- ifelse(hh$nblayoffs>5,simplelag(hh$playoff,by=hh$est,outside=0),hh$playoff)
hh$flayoff <- ifelse(hh$nblayoffs>5,forward(hh$flayoff,by=hh$est,outside=0),hh$flayoff)

# hh$playoff_sz <- retain(hh$nblayoffs,hh$nblayoffs>0 | hh$est!=simplelag(hh$est),outside=0)
# hh$flayoff_sz <- obtain(hh$nblayoffs,hh$nblayoffs>0  | hh$est != forward(hh$est),outside=0)


# Stacking data together
gc()
fff <- NULL
for (y in c(2002:2014)) {
  hh$treated_y <- ( hh$year<y+5 &  hh$year>y-5 &
                      ((hh$nblayoffs>5)*(hh$year==y) | hh$playoff==y | hh$flayoff==y) &
                      (hh$flayoff==0 | ( ( (hh$flayoff>y)*hh$year ) < hh$flayoff )) &
                      (hh$playoff==0 | (hh$year > ((hh$playoff<y)*hh$playoff)))  &
                      ((hh$nblayoffs>5)*(hh$year != y) ==0)
  )*1
  hh$treated_y_sum <- ave(hh$treated_y,hh$est,FUN=sum)
  
  hh$control_y <- ( hh$treated_y_sum==0 &
                      hh$year<y+5 &  hh$year>y-5 &
                      (hh$nblayoffs<6) &
                      (hh$flayoff==0 | (hh$year < hh$flayoff) ) &
                      (hh$playoff==0 | (hh$year > hh$playoff  )) 
  )*1
  
  ff <- hh[hh$treated_y %in% 1 | hh$control_y %in% 1,
           c("est","year","count","lnnbwkrs_cumneg","f9010xf9010","f9010_w","f9010_w2","nblayoffs",
             "playoff","flayoff","r_lyof","treated_y","control_y")]
  ff$eyear <- y
  ff$y_nblayoffs <- ave((ff$nblayoffs>0)*ff$year,ff$est,FUN=function(x) max(x,na.rm=T))
  ff$min_playoff <- ave(ff$playoff,ff$est,FUN=function(x) min(x,na.rm=T))
  ff$nb_obs <- ave(ff$year>0,ff$est,FUN=sum)
  ff$min_year <- ave(ff$year,ff$est,FUN=min)
  ff$max_year <- ave(ff$year,ff$est,FUN=max)
  fff <- rbind(fff,ff[ff$nb_obs>1 &  ff$min_year<=y & ff$max_year>=y,])
}



table(fff$y_nblayoffs,fff$treated_y)
nrow(fff)

table(fff$nblayoffs>0,fff$year)

#-----------------------------------------------------------------------------#
#### MODELS 5+ layoffs ####
#-----------------------------------------------------------------------------#
gc()
fff$firm <- substr(fff$est,1,9)

mod90_3 <- felm(I(100*f9010xf9010) ~ factor(paste(eyear,year))
                + I((year - y_nblayoffs) %in% -4 )
                + I((year - y_nblayoffs) %in% -3 )
                + I((year - y_nblayoffs) %in% -2 )
                # + I((year - y_nblayoffs) %in% -1 & y_nblayoffs==eyear)
                + I((year - y_nblayoffs) %in% 0 )
                + I((year - y_nblayoffs) %in% 1 )
                + I((year - y_nblayoffs) %in% 2 )
                + I((year - y_nblayoffs) %in% 3 )
                + I((year - y_nblayoffs) %in% 4 )
                |paste(eyear,est)|0|firm, 
                weights=fff$f9010_w2, 
                data=fff  )
summary(mod90_3)
screenreg(list(mod90_3),stars=c(0.1,0.05,0.01),digits=3)
htmlreg(list(mod90_3),stars=c(0.1,0.05,0.01),file="did_layoff5_stacked.html",digits=3)



mod90_3 <- felm(I(100*f9010xf9010) ~ factor(paste(eyear,year))
                + I((year - y_nblayoffs) %in% -4 )
                + I((year - y_nblayoffs) %in% -3 )
                + I((year - y_nblayoffs) %in% -2 )
                # + I((year - y_nblayoffs) %in% -1 & y_nblayoffs==eyear)
                + I((year - y_nblayoffs) %in% 0 )
                + I((year - y_nblayoffs) %in% 1 )
                + I((year - y_nblayoffs) %in% 2 )
                + I((year - y_nblayoffs) %in% 3 )
                + I((year - y_nblayoffs) %in% 4 )
                + log(count)
                + lnnbwkrs_cumneg
                |paste(eyear,est)|0|firm, 
                weights=fff$f9010_w2, 
                data=fff  )
summary(mod90_3)
screenreg(list(mod90_3),stars=c(0.1,0.05,0.01),digits=3)
htmlreg(list(mod90_3),stars=c(0.1,0.05,0.01),file="did_layoff5_stacked_size.html",digits=3)


